C.................... AURSUB.FOR ................................. 
C--- This subroutine reads in auroral electron production frequencies
C--- and converts them to production rates for the FLIP model. Since
C--- the auroral spatial grid differs from both the field line grid
C--- and the vertical grid used by the FLIP model, interpolation is
C--- required. This program was written by P. Richards in August 1991
      SUBROUTINE AURSUB(MINJZ,MAXJZ,TIME,ZALT,ON,O2N,N2N,N,TI
     > ,EHEAT,DTMAX,ZWR)
      DOUBLE PRECISION EHEAT(3,401),ZALT(401),ON(401),O2N(401),N2N(401)
     >  , N(4,401),TI(3,401),AURALT,DTMAX
      INTEGER FEXISTS,JCALL
      REAL XN(3,401),XNE(401),XTE(401)
      REAL ISAV,ESAV
      REAL TIME
      !.. two commons to hold the EUV and AURoral production rates
      COMMON/EUVPRD/EUVION(3,12,401),PEXCIT(3,12,401),PEPION(3,12,401)
      COMMON/AURPRD/AURION(3,12,401),AUREX(3,12,401),AURALT(401),IAMAX
      COMMON/AURP2S/ISAV(3,12,401),ESAV(3,12,401),AURZ(401),IAMX
      COMMON/BPRAUR/PAUION(3,12,401),PAUEXC(3,12,401)
      COMMON/IMA/IM,IA

      !.. MAXIND is used for storing in conjugate portion of the arrays. It
      !.. assumes that MAXJZ=JMAX/2 where JMAX = # of field line grid points
      MAXIND=2*MAXJZ+1

      DO 2 IS=1,3
      DO 2 IK=1,12
      DO 2 IJ=MINJZ,MAXIND
        AURION(IS,IK,IJ)=0.0
        AUREX(IS,IK,IJ)=0.0
        PAUION(IS,IK,IJ)=0.0
        PAUEXC(IS,IK,IJ)=0.0
  2   CONTINUE

      JCALL=JCALL+1  !.. count number of times routine is called
      !... see if auroral file exists
      IF(JCALL.EQ.1) CALL TFILE(14,FEXISTS)
      IF(FEXISTS.EQ.0) RETURN

      CALL AURORA(TIME,ZWR,AURFAC,DTMAX_S)

      DTMAX=DTMAX_S  !.. set maximum time step for aurora

      IF(AURFAC.LE.0.001) RETURN

      !... Load FLIP densities into XN to aid looping
      DO 5 J=MINJZ,MAXIND
        XN(1,J)=ON(J)
        XN(2,J)=O2N(J)
        XN(3,J)=N2N(J)
        XTE(J)=TI(3,J)
        XNE(J)=N(1,J)+N(2,J)+N(3,J)+N(4,J)
        !.. WRITE(116,145) J,ZALT(J),EHEAT(3,J)
        !.. >    ,EUVION(1,1,J),PEXCIT(1,1,J),PEPION(1,1,J)
 5    CONTINUE
      !--- adjust max alt index down to avoid divide by zero.
      DO 6 J=MINJZ,MAXJZ
        IF(ZALT(J).LT.1500) JALT=J
 6    CONTINUE
      MAXJZ=JALT-1

      !... transfer real altitude to double altitude for 
      !... interpolation routines
      IAMAX=IAMX
      DO J=1,IAMAX
        AURALT(J)=AURZ(J)
      ENDDO

      !-- modify input frequencies by AURFAC
      DO 9 IK=1,12
      DO 9 IS=1,3
      DO 9 J=1,IAMAX
         AURION(IS,IK,J)= AURFAC * ISAV(IS,IK,J)
         AUREX(IS,IK,J)= AURFAC * ESAV(IS,IK,J)
 9    CONTINUE

      !.. Interpolate auroral production to FLIP field line grid 
      !.. The same frequencies are used for the conjugate point using ICONJ
      !.. However, they are multiplied by the appropriate density
      DO 220 ICONJ=1,2
         DO 10 IK=1,12
         DO 10 IS=1,3
         DO 10 JJ=MINJZ,MAXJZ
            J=JJ
            IF(ICONJ.EQ.2) J=MAXIND+1-JJ
            !... make sure FLIP altitude is inside auroral grid
            IF(AURALT(1).LE.ZALT(J).AND.AURALT(IAMAX).GE.ZALT(J)) THEN
              PAUION(IS,IK,J)=
     >          FTEUVP(IM,IS,IK,1,IAMAX,AURALT,AURION,ZALT(J))*XN(IS,J)
              PAUEXC(IS,IK,J) =
     >          FTEUVP(IM,IS,IK,1,IAMAX,AURALT,AUREX,ZALT(J))*XN(IS,J)
            ENDIF
 10      CONTINUE

      !.. electron heating and miscellaneous production rates
         DO 15 II=MINJZ,MAXJZ
            I=II
            IF(ICONJ.EQ.2) I=MAXIND+1-II
            EHEAT(3,I)=EHEAT(3,I)+PAUEXC(1,9,I)*(XNE(I)**0.97)/XN(1,I)
            PAUEXC(3,11,I)=PAUEXC(3,11,I)/XN(3,I)
      !..  add in 1356 from O2 dissociation
            PAUEXC(1,6,I)= PAUEXC(1,6,I)  + PAUEXC(2,9,I)
 15      CONTINUE
C
          !-----debug writes
          IF(IAMAX.LT.0) GO TO 220
          IF(AURFAC.EQ.-9999) THEN
           DO 324 J=MINJZ,MAXJZ
             WRITE(116,145) J, ZALT(J), EHEAT(3,J),PAUEXC(1,3,J)
     >         ,PAUEXC(2,3,J),PAUEXC(3,3,J),PAUION(1,1,J)
     >         ,PAUION(2,1,J),PAUION(3,1,J)
 324       CONTINUE
           DO 325 J=MINJZ,MAXJZ
              WRITE(116,145) J, ZALT(J),(PAUEXC(1,K,J),K=1,9)
     >           ,(PAUION(1,K,J),K=1,3)
 325       CONTINUE
           DO 326 J=MINJZ,MAXJZ
              WRITE(116,145) J, ZALT(J),(PAUEXC(2,K,J),K=1,9)
     >           ,(PAUION(2,K,J),K=1,3)
 326       CONTINUE
           DO 327 J=MINJZ,MAXJZ
              WRITE(116,145) J, ZALT(J),(PAUEXC(3,K,J),K=1,9)
     >           ,(PAUION(3,K,J),K=1,3)
 327       CONTINUE
 145       FORMAT(1X,I3,F7.1,1P,22E9.1)
          ENDIF

 220  CONTINUE

 333    RETURN

      !.. error trap
 298  CONTINUE
      WRITE(6,*) '   *** ERROR: in Aurora (unit FOR014) data file'
      CALL RUN_ERROR    !.. print ERROR warning in output files
      STOP
      END
C:::::::::::::::::::::::: BLOCK DATA AURBLK :::::::::::::::
      BLOCK DATA AURBLK
      COMMON/IMA/IM,IA
      DATA  IM/0/, IA/0/
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE PRODIT(J,ISW,MINJ,MAXJ)
C...... this program interpolates the various production rates for
C...... the FLIP model vertical grid and saves them so that they
C...... do not have to be recalculated
C...... This program was written by P. Richards in August 1991
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL ENTION,PNTXIT,PNTION,ANTION,ANTEX
      DIMENSION  VIBFAC(7)
      COMMON/ENLOSS/ENL(20,VD)
      COMMON/NPRD/COLUM(3,401),OTHPR1(6,401),OTHPR2(6,401),OTHPR3(6,401)
      COMMON/PXS/VK(6),PXN2A(VD),PEEX(7,VD),PRN2P(7,VD),JVMAX
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/EPRDNT/ENTION(3,12,VD),PNTXIT(3,12,VD),PNTION(3,12,VD)
      COMMON/APRDNT/ANTION(3,12,VD),ANTEX(3,12,VD)
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/PRTER/PUVN2P(VD),PUVO2D(VD),PEO1D(VD),PUVO2P(VD)
     >  ,PUOP2P(VD),UVDISN(VD),PSRNO(VD),PO1DSR(VD),PO1SEL(VD)
     >  ,PN2A3S(VD),PDISOP(VD),PUVIOP(VD),PRN2XS(VD),PRN2AP(VD)
     >  ,PRN2BS(VD),DISNP(VD),DISN4S(VD),DISN2D(VD),DISN2P(VD),PRHEP(VD)
     >  ,HEPLUS(VD),NPLUS(VD),PLYNOP(VD),PHOTN(VD),HE10830(VD)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)
      COMMON/FJS/N(4,401),TI(3,401),FF(20)
      DATA VIBFAC/1.,.36,.2,.13,.11,.21,.0/
      COMMON/IMA/IM,IA

      !.. production of n2+(x2s),n2+(a2p),n2+(b2s) by EUV, e*, and aurora
      PRN2XS(J)= ENTION(3,1,J)+PNTION(3,1,J)+ANTION(3,1,J)
      PRN2AP(J)= ENTION(3,2,J)+PNTION(3,2,J)+ANTION(3,2,J)
      PRN2BS(J)= ENTION(3,3,J)+PNTION(3,3,J)+ANTION(3,3,J)
      PUVN2P(J)=PRN2XS(J)+PRN2AP(J)+PRN2BS(J)
      !.. production of O+ states - includes contribution from 2P*, 4P
      PUVIOP(J)= ENTION(1,7,J)+PNTION(1,7,J)+ANTION(1,7,J)
      PUVO2D(J)= ENTION(1,8,J)+PNTION(1,8,J)+ANTION(1,8,J)
      PUOP2P(J)= ENTION(1,9,J)+PNTION(1,9,J)+ANTION(1,9,J)
      !.. production of O(1D) and O(1S) by e* and aurora
      PEO1D(J)= PNTXIT(1,1,J)+ANTEX(1,1,J)
      PO1SEL(J)= PNTXIT(1,2,J)+ANTEX(1,2,J)
      !.. total production of O2+  - includes all non-diss states
      PUVO2P(J)= ENTION(2,7,J)+PNTION(2,7,J)+ANTION(2,7,J)
      !... dissociation of N2 to yield N+
      DISNP(J)= ENTION(3,4,J)+ENTION(3,5,J)+ENTION(3,6,J)
     > +PNTION(3,6,J)+ANTION(3,6,J)
      !.. Total dissociation of N2 to give metastable nitrogen atoms
      DISN4S(J)=PNTXIT(3,5,J)+ANTEX(3,5,J)+DISNP(J)/2.0
      DISN2D(J)=PNTXIT(3,6,J)+ANTEX(3,6,J)+DISNP(J)/2.0
      DISN2P(J)=PNTXIT(3,7,J)+ANTEX(3,7,J)
      !.... photodissociation of n2 into n(2d) and n(4s) and n(2p)
      UVDISN(J)=FTX(1,MINJ,MAXJ,Z,OTHPR1,ZR(J))
      !.. production of He+
      PRHEP(J)=FTX(2,MINJ,MAXJ,Z,OTHPR1,ZR(J))
      ENL(19,J)=FTX(4,MINJ,MAXJ,Z,OTHPR1,ZR(J))
      !... production of o(1d) by schumann-runge
      PO1DSR(J)=FTX(3,MINJ,MAXJ,Z,OTHPR1,ZR(J))
      ENL(11,J)=FTX(5,MINJ,MAXJ,Z,OTHPR1,ZR(J))
      !.... total production of N2(A) by electron impact excitation
      PN2A3S(J)= PNTXIT(3,1,J) + PNTXIT(3,2,J)+ PNTXIT(3,3,J)
     > + ANTEX(3,1,J) + ANTEX(3,2,J)+ ANTEX(3,3,J)
      !.... production frequency for e* excitation of N2(A)
      PXN2A(J)=PN2A3S(J)/N2N(J)
      !.. dissociative ionization of o2 - fix e* and auroral
      PDISOP(J)= ENTION(2,4,J)+ENTION(2,5,J)+PNTION(2,4,J)+ANTION(2,4,J)
     > +PNTION(2,5,J)+ANTION(2,5,J)

      !.. direct excit by e* (<5eV) using small x-s.
      PHEX= (PNTXIT(3,8,J) + ANTEX(3,8,J))/N2N(J)
      DO 32 IE=1,7
 32   PEEX(IE,J)=PHEX*VIBFAC(IE)
      !.. He+ and N+
      HEPLUS(J)=FITDEN(4,MINJ,MAXJ,Z,N,ZR(J))
      NPLUS(J)=FTXION(1,MINJ,MAXJ,Z,XIONN,ZR(J))
      !..-- hv + NO -> NO+ + e
      PLYNOP(J)=FTX(1,MINJ,MAXJ,Z,OTHPR2,ZR(J))
      !..-- hv (SRB) + NO -> N(4S) + O
      PSRNO(J)=FTX(2,MINJ,MAXJ,Z,OTHPR2,ZR(J))
      !..-- hv + N(4S) -> N+ + e
      PHOTN(J)=FTX(3,MINJ,MAXJ,Z,OTHPR2,ZR(J))
       !..-- e* + He -> He 10830 (No aurora at this time)
      He10830(J)=FTX(4,MINJ,MAXJ,Z,OTHPR2,ZR(J))
      RETURN
      END
C::::::::::::::::::::::::: VNTERP ::::::::::::::::::::::::::::::::::::::
C
C...SUBROUTINE: VNTERP, MODULE:     AURSUB.FOR
C...PROGRAMMER: G. GERMANY, DATE:       11/4/91
C...PURPOSE:    TO INTERPOLATE PRODUCTION RATES FROM THE FIELD
C...            LINE AND AURORAL ALTITUDE GRIDS TO THE VERTICAL
C...            ALTITUDE GRID
C...NOTES:      PRODUCTION RATES ARE PASSED IN BY COMMON BLOCKS
C...            EUVPRD AND AURPRD.  THE INTERPOLATED VALUES ARE
C...            RETURN IN COMMON BLOCKS EPRDNT AND APRDNT
C
C...            THE ROUTINE IS CALLED FROM A DO LOOP WITH COUNTER
C...            JC.  IT RUNS ONLY ONCE IN THE DO LOOP.  IF THE
C...            COUNTER JC>1 THEN ROUTINE RETURNS
C...VARIABLES:  SEE BELOW
C... DEBUG  --> LOGICAL FOR DEBUGGING
C... IUN    --> UNIT NUMBER FOR DEBUG OUTPUT
C... AURALT --> AURORAL         ALTITUDE ARRAY /IAMAX PTS/
C... ZR     --> FLIP VERTICAL   ALTITUDE ARRAY /  VD  PTS/
C... Z      --> FLIP FIELD LINE ALTITUDE ARRAY /  60  PTS/
C...            PE IS SHORTHAND FOR 'PHOTO-ELECTRON'
C... EUVION --> PRODUCTION RATE ARRAY DUE TO EUV IONIZATION
C... PEXCIT --> PRODUCTION RATE ARRAY DUE TO PE  EXCITATION
C... PEPION --> PRODUCTION RATE ARRAY DUE TO PE  IONIZATION
C... ENTION --> EUVION INTERPOLATED TO VERTICAL ALTITUDE GRID
C... PNTXIT --> PEXCIT INTERPOLATED TO VERTICAL ALTITUDE GRID
C... PNTION --> PEPION INTERPOLATED TO VERTICAL ALTITUDE GRID
C... AURION --> PRODUCTION RATE ARRAY DUE TO AURORAL IONIZATION
C... AUREX  --> PRODUCTION RATE ARRAY DUE TO AURORAL EXCITATION
C... ANTION --> AURION INTERPOLATED TO VERTICAL ALTITUDE GRID
C... ANTEX  --> AUREX  INTERPOLATED TO VERTICAL ALTITUDE GRID
C-----------------------------------------------------------------------
      SUBROUTINE VNTERP(JB,MINJ,MAXJ)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL EUVION,PEXCIT,PEPION,AURION,AUREX,ENTION,PNTXIT,PNTION
     &  ,ANTION,ANTEX,TOTION,TOTEX
      real FTEUVP            !!!! f90 added
      LOGICAL*4 DEBUG
      DATA DEBUG/.FALSE./IUN/50/

      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      !.. original production rates
      COMMON/EUVPRD/EUVION(3,12,401),PEXCIT(3,12,401),PEPION(3,12,401)
      COMMON/AURPRD/AURION(3,12,401),AUREX(3,12,401),AURALT(401),IAMAX
      !.. interpolated production rates
      COMMON/EPRDNT/ENTION(3,12,VD),PNTXIT(3,12,VD),PNTION(3,12,VD)
      COMMON/APRDNT/ANTION(3,12,VD),ANTEX(3,12,VD)
      COMMON/PTOTV/TOTION(3,12,VD),TOTEX(3,12,VD)

      !..- interpolate production rates and frequencies to vertical grid
      DO 10 IS=1,3
      DO 10 IK=1,12
      DO 10 J=1,JB+1
         ENTION(IS,IK,J)=FTEUVP(0,IS,IK,MINJ,MAXJ,Z,EUVION,ZR(J))
         PNTXIT(IS,IK,J)=FTEUVP(0,IS,IK,MINJ,MAXJ,Z,PEXCIT,ZR(J))
         PNTION(IS,IK,J)=FTEUVP(0,IS,IK,MINJ,MAXJ,Z,PEPION,ZR(J))
         ANTION(IS,IK,J)=FTEUVP(0,IS,IK,1,IAMAX,AURALT,AURION,ZR(J))
         ANTEX(IS,IK,J) =FTEUVP(0,IS,IK,1,IAMAX,AURALT,AUREX,ZR(J))
 10    CONTINUE

      !.. Multiply the auroral production frequencies by density
      DO 20 IS=1,3
      DO 20 J=1,JB+1
         IF(IS.EQ.1) DENS = ON(J)
         IF(IS.EQ.2) DENS = O2N(J)
         IF(IS.EQ.3) DENS = N2N(J)
         DO 20 IK=1,12
            ANTION(IS,IK,J)=ANTION(IS,IK,J)*DENS
            ANTEX(IS,IK,J) =ANTEX(IS,IK,J)*DENS
 20   CONTINUE

      DO 30 J=1,JB+1
      !..     electron heating rate take [O] out again
         ANTEX(1,9,J)=ANTEX(1,9,J)/ON(J)
      !..    add O2 dissociation of 1356
         ANTEX(1,6,J)=ANTEX(1,6,J)+ANTEX(2,9,J)
         !-- extract density from total auroral input and 1200 prod because
         !-- they already had densities associated with them
         ANTION(1,11,J)=ANTION(1,11,J)/ON(J)
         ANTEX(3,11,J)=ANTEX(3,11,J)/N2N(J)
      !..    add contributions from 2 highest O+ metastables to 3 lowest
         ENTION(1,7,J) = ENTION(1,1,J) + ENTION(1,4,J)
         ENTION(1,8,J) = ENTION(1,2,J) + ENTION(1,5,J)/1.3
         ENTION(1,9,J) = ENTION(1,3,J) + ENTION(1,5,J)/4.3

         PNTION(1,7,J) = PNTION(1,1,J) + PNTION(1,4,J)
         PNTION(1,8,J) = PNTION(1,2,J) + PNTION(1,5,J)/1.3
         PNTION(1,9,J) = PNTION(1,3,J) + PNTION(1,5,J)/4.3

         ANTION(1,7,J) = ANTION(1,1,J) + ANTION(1,4,J)
         ANTION(1,8,J) = ANTION(1,2,J) + ANTION(1,5,J)/1.3
         ANTION(1,9,J) = ANTION(1,3,J) + ANTION(1,5,J)/4.3
 
      !.. SUM non-diss. states to get total O2+ and N2+ production
         ENTION(2,7,J)= ENTION(2,1,J)+ENTION(2,2,J)+ENTION(2,3,J)
         PNTION(2,7,J)= PNTION(2,1,J)+PNTION(2,2,J)+PNTION(2,3,J)
         ANTION(2,7,J)= ANTION(2,1,J)+ANTION(2,2,J)+ANTION(2,3,J)

         ENTION(3,7,J)= ENTION(3,1,J)+ENTION(3,2,J)+ENTION(3,3,J)
         PNTION(3,7,J)= PNTION(3,1,J)+PNTION(3,2,J)+PNTION(3,3,J)
         ANTION(3,7,J)= ANTION(3,1,J)+ANTION(3,2,J)+ANTION(3,3,J)
 30   CONTINUE

      !.. Total ionization and excitation rates
      DO 40 IS=1,3
      DO 40 IK=1,12
      DO 40 J=1,JB+1
         TOTION(IS,IK,J) = ENTION(IS,IK,J)+PNTION(IS,IK,J)
     >       +ANTION(IS,IK,J)
         TOTEX(IS,IK,J) = PNTXIT(IS,IK,J)+ANTEX(IS,IK,J)
 40   CONTINUE

      RETURN
      END
C:::::::::::::::::::::::::: FTXION :::::::::::::::::::::::::::::::::
C...... ftn and ftt are interpolating programs
      FUNCTION FTXION(I,JMIN0,JMAX0,ZP,T,Z)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION ZP(401),T(9,401)

      !.. Use bisection to find the nearest altitude ZP(K) to desired
      !.. altitude (Z)
      CALL BISPLT(JMIN0,JMAX0,Z,ZP,I1)

      !.. Now interpolate to get desired production rate. Test to make
      !.. sure we haven't run off the end of the field line grid
 2    IF(I1.EQ.JMAX0) THEN
         FTXION=T(I,I1)
      ELSE IF(T(I,I1).LT.0.OR.T(I,I1+1).LT.0) THEN
         FTXION=0.0
      ELSE
         FACTOR=(Z-ZP(I1))/(ZP(I1+1)-ZP(I1))
         FTXION=T(I,I1)+FACTOR*(T(I,I1+1)-T(I,I1))
      ENDIF
 4    RETURN
      END
C::::::::::::::::::::::::::: FTX ::::::::::::::::::::::::::::::::::
C.......exponential interpolation of production rates (fitting y= A exp(Bx)
      FUNCTION FTX(I,JMIN0,JMAX0,ZP,T,Z)
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      DIMENSION ZP(401),T(6,401)

      !.. Use bisection to find the nearest altitude ZP(K) to desired 
      !.. altitude (Z)
      CALL BISPLT(JMIN0,JMAX0,Z,ZP,K)

      !.. Now interpolate to get desired production rate. Test to make
      !.. sure we haven't run off the end of the field line grid
 2    IF(K.EQ.JMAX0) THEN
         FTX=T(I,K)
      ELSE IF(T(I,K).LT.0.OR.T(I,K+1).LT.0) THEN
         FTX=0.0
      ELSE 
         TK=(T(I,K)) + 1.0E-22
         TK1=(T(I,K+1)) + 1.0E-22
         B=DLOG(TK/TK1)/(ZP(K)-ZP(K+1))
         FTX=TK*DEXP(B*(Z-ZP(K)))
         IF(K.EQ.1.AND.Z.EQ.ZP(K)) FTX=TK
      ENDIF
      RETURN
      END
C:::::::::::::::::::::::::: FITDEN :::::::::::::::::::::::::::::::::
C..exponential interpolation of densities (fitting y= A exp(Bx)
      FUNCTION FITDEN(I,JMIN0,JMAX0,ZP,T,Z)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION ZP(401),T(4,401)

      !.. Use bisection to find the nearest altitude ZP(K) to desired 
      !.. altitude (Z)
      CALL BISPLT(JMIN0,JMAX0,Z,ZP,K)

      !.. Now interpolate to get desired production rate. Test to make
      !.. sure we haven't run off the end of the field line grid
 2    IF(K.EQ.JMAX0) THEN
         FITDEN=T(I,K)
      ELSE IF(T(I,K).LT.0.OR.T(I,K+1).LT.0) THEN
         FITDEN=0.0
      ELSE 
         TK=(T(I,K)) + 1.0E-22
         TK1=(T(I,K+1)) + 1.0E-22
         B=DLOG(TK/TK1)/(ZP(K)-ZP(K+1))
         FITDEN=TK*DEXP(B*(Z-ZP(K)))
      ENDIF
      RETURN
      END
C:::::::::::::::::::::::: FITTEM ::::::::::::::::::::::::::::::::::::
C...... Linear interpolation for temperatures (y = A.x + B)
      FUNCTION FITTEM(I,JMIN0,JMAX0,ZP,T,Z)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION ZP(401),T(3,401)

      !.. Use bisection to find the nearest altitude ZP(K) to desired 
      !.. altitude (Z)
      CALL BISPLT(JMIN0,JMAX0,Z,ZP,I1)

      !.. now do linear fit
 2    IF(I1.EQ.JMAX0) THEN
         FITTEM=T(I,I1)
      ELSE IF(T(I,I1).LT.0.OR.T(I,I1+1).LT.0) THEN
         FITTEM=0.0
      ELSE
         FACTOR=(Z-ZP(I1))/(ZP(I1+1)-ZP(I1))
         FITTEM=T(I,I1)+FACTOR*(T(I,I1+1)-T(I,I1))
      ENDIF
 4     RETURN
       END
C::::::::::::::::::::::::::: FTEUVP ::::::::::::::::::::::::::::::::::
C...... exponential interpolation for EUV production rates.
C...... This program was written by P. Richards in August 1991
      FUNCTION FTEUVP(IM,I,J,JMIN0,JMAX0,ZP,T,Z)
      IMPLICIT REAL(A-H,L,N-Z)
      DOUBLE PRECISION ZP,Z
      DIMENSION ZP(401),T(3,12,401)
      !.. Check to ensure Z is in the range of altitudes
      FTEUVP=0.0
      IF(JMAX0.LE.0) RETURN

      !.. Use bisection to find the nearest altitude ZP(K) to desired 
      !.. altitude (Z)
      CALL BISPLT(JMIN0,JMAX0,Z,ZP,K)

      !.. added to stop running past the upper boundary 2008-12-06
      IF(K+1.GT.JMAX0) K=K-1
      !.. Now interpolate to get desired production rate
 2    IF(K.EQ.JMAX0) THEN
         FTEUVP=T(I,J,K)
      ELSE IF(T(I,J,K).LT.0.OR.T(I,J,K+1).LT.0) THEN
         FTEUVP=0.0
      ELSE 
         TK=(T(I,J,K)) + 1.0E-22
         TK1=(T(I,J,K+1)) + 1.0E-22
         B=ALOG(TK/TK1)/(ZP(K)-ZP(K+1))
         FTEUVP=TK*EXP(B*(Z-ZP(K)))
         IF(K.EQ.1.AND.Z.EQ.ZP(K)) FTEUVP=TK
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::: BISPLT ::::::::::::::::::::::::::::::::::
C--- Use bisection to find the nearest altitude ZP(K) to desired
C--- altitude (Z). Modified in Feb 93 to do both the northern and
C--- southern hemispheres.
      SUBROUTINE BISPLT(BEGIN,JMAX0,Z,ZP,FIRST)
      IMPLICIT REAL(A-H,L,N-Z)
      DOUBLE PRECISION ZP(401),Z
      INTEGER FIRST,LAST,MID,BEGIN

         FIRST=BEGIN
         LAST=JMAX0
         ITER=0
      !.. This section for northern hemisphere
      IF(ZP(BEGIN).LE.ZP(JMAX0)) THEN
 10      CONTINUE
             MID=(FIRST+LAST)/2
             IF(Z.LE.ZP(MID)) THEN
                LAST=MID-1
             ELSE
                FIRST=MID+1
             ENDIF  
         IF(FIRST.LT.LAST) GO TO 10
         IF(Z.LT.ZP(FIRST)) FIRST=FIRST-1
         IF(Z.LT.ZP(BEGIN)) FIRST=BEGIN

      !.. Southern hemisphere is different because altitude decreasing
      ELSE
 110      CONTINUE
             MID=(FIRST+LAST)/2
             IF(Z.GE.ZP(MID)) THEN
                LAST=MID-1
             ELSE
                FIRST=MID+1
             ENDIF  
             ITER=ITER+1
         IF(FIRST.LT.LAST) GO TO 110
         IF(Z.GT.ZP(FIRST)) FIRST=FIRST-1
         IF(Z.LT.ZP(JMAX0)) FIRST=JMAX0
      ENDIF

 901  FORMAT(1X,4I4,1P,9E10.2)
      RETURN
      END
